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1. Introduction 

In the experimental exploration of complex systems, such es those encountered in 
life-science, geology, or astronomy, it is not unusual that the experimenter discovers 
oscillations of unknown origin in a measured time-series. The experimenter would then 
usually try to characterize these oscillations in a form that admits an identification 
of their source - the oscillator. The conventional first choice is a characterization of 
the oscillations by their "frequency". For an ideal, periodically oscillating signal x(t), 
the smallest number T > such that x(t) — x(t + T) for all t is the period of the 
signal, and its (angular) frequency is defined by ui — 2ir/T. There is a good reason 
for choosing this particular characterization. All other properties of an ideal, periodic 
signal, i.e. its waveform and amplitude, are subject to distortions along the signal 
pathway from the oscillator to the detector. In fact, linear filtering along the signal 
pathway would generally be sufficient to modify the waveform and the amplitude in 
an arbitrary way. And the precise properties of the signal pathway are unknown in the 
setting considered here. The frequency information is the only sure fact. For ideal, 
periodic oscillations, these observations are too obvious to deserve much discussion. 
But for non-ideal oscillations, as they are frequently encountered in complex systems, 
the situation is less clear. 

A large variety of methods is being used to determine a "frequency" for non-ideal 
oscillations, and not all of them are equally robust to filtering and other distortions. 
Two major kinds of methods can be distinguished: Firstly, there are period- counting 
methods, where, from the number of oscillation periods n(At) in a time interval 
[to , to + At] , the frequency is determined as 

r 27m ( At ) m 

^count = hm — . (1) 

At— >oo LAX 

(Finite sample-size effects are not discussed here.) The methods differ in the criteria 
used for counting individual periods (e.g., local maxima, zero-crossings). Secondly, 
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there are spectral methods, where the frequency Lo spcc characterizes the position of 
a peak in an estimate of the power spectral density S x (v) P [2] of the signal x{t) 
(e.g., EJ 00 ED EI)- Often, the frequency with maximum power w pca k is used. 
The term spectral methods shall here also include methods based on estimates of 
the autocorrelation function of the signal, since this is related to the spectral density 
by a simple Fourier transformation. 

For weak distortions and not too irregular oscillations, period-counting methods 
can be just as unequivocal as frequency measurements for ideal oscillations. This 
is why they are routinely used in high-precision frequency (or time) measurements. 
They are also naturally associated with mode- locking phenomena [S]. But for stronger 
distortions and more irregular oscillations, this robustness is reverted to its contrary. 
Unequivocally identifying individual periods of oscillation then becomes difficult. In 
these situations spectral methods are generally preferred. However, it is obvious that 
spectral methods are not robust to filtering along the signal pathway either. By linear 
filtering, the power spectrum can be modified nearly arbitrarily. 

How much can the concept of period-counting frequency measurement be 
extended to distorted time series? A partial answer is given in this work. In section l2~D 
a generalized period-counting frequency measure, the topological frequency, is defined. 
It is based on the approximate reconstruction of the phase-space trajectory of the 
oscillator. In section 12.21 it is shown that the topological frequency is robust with 
respect to nearly arbitrary linear filtering. This has three important implications: (i) 
At least as long as the signal pathway acts as a linear filter, the topological frequency 
is a characteristic of the (typically nonlinear) oscillator alone, (ii) Filtering of the 
signal, in order to remove noise and other undesirable components, does no harm to 
the result for the frequency. In view of (i) and (ii), the topological frequency can be 
considered to be robust with respect to both kinds of distortions, filtering and noise. 
Finally, (iii) the results of frequency measurements using spectral methods can deviate 
arbitrarily from the topological frequency. This point is made rigorous in section \2. 31 

Not for all oscillatory time series can the topological frequency be defined. In 
particular, linear time series driven by Gaussian noise are excluded. Weakly nonlinear 
models for noisy time series can interpolate between linear Gaussian oscillations and 
ideal periodicity. For these models, the influence of filtering on a weaker period- 
counting frequency measure, the average or phase frequency (see below), is investigated 
in sectional The susceptibility of the phase frequency to filtering is found to decay 
rapidly with the degree of nonlincarity. Section 01 contains some concluding remarks. 

2. Topological Frequency 

2.1. Definitions 

The theory becomes more transparent in a discrete-time representation. Let {xt} be 
an infinite, real-valued time series sampled at equally spaced times starting at t = 0. 
Measure time in units of the sampling interval. Define the spectral density S x (lo) of 
{x t } as 



where (-) t denotes temporal averaging (t > 0). 

Let the trajectory p(t) of a time series {x{\ in A-dimensional delay space be 
defined by p(t) — (xt, %t-ii ■ ■ ■ > &t—N+i) for integer t and by linear interpolation^: for 

\ In practice, higher order interpolation might sometimes be useful. 




(2) 



T — — 



OO 
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non-integer t. Frequency will here be defined with respect to a Poincare section or 
counter, which is an (TV — l)-dimcnsional, oriented manifold M with boundary dM 
and interior Int M = M\dM, embedded in the TV-dimensional delay space. 

Let n{t\) be the oriented number of transitions of the trajectory p{t) through 
Int M in the time interval < t < t\. That is, a transition through IntM in positive 
(negative) direction increments (decrements) n{t\) by one. For example, a positive- 
slope zero-crossing counter in 2-dimensional delay space would be given by§ 

M = {{v u v 2 ) e M 2 \v 1 +v 2 = 0, V! > v 2 }. (3) 

Define the topological frequency uj m _ x of {x t } with respect to a counter M, as 

lo m ,x := hm — (4) 

t^QO t 

provided the limit exists and there is a d > such that p(t) has for alH > a distance 
> d from dM. By construction, lom,x is invariant under not too large deformations 
of {xt} and M. For example, if the trajectory contains a loop which comes close 
to IntM but does not encircle dM, a small deformation of {x t } or M might make 
the loop intersect IntM. But, since this intersection comprises two transitions of the 
trajectory through Int M, one of which is positive and one of which is negative, the 
value of n(t) does not change for large enough t. Configurations with p(t) tangential to 
Int M can be evaluated as either of both limiting cases - intersecting or not - without 
effecting ojm,x- Drastically different counters can lead to different frequencies. But 
each of these is sharply defined. 

2.2. Invariance under filtering 

It can be shown that for any bounded time series {a;*} the topological frequency is 
invariant under nearly arbitrary linear filtering: 

Theorem 1. Let {y t } be obtained from a bounded time series {x t } by linear, causal 
filtering, 

oo 

y t := ^a k x t -k- (5) 

fc=0 

Assume that, for some r > I, 



< 



^a k z k 



fe=0 



< oo (6) 



for all complex z, r^ 1 < \z\ < r. (This excludes, for example, filters which fully 
block some frequencies.) Let M be a counter and ujm,x be defined. Then there is, at 
sufficiently high embedding dimension, a counter M' such that ojm', v is defined and 

UM',y — Um,x- 

Proof. This is most easily seen by the following explicit construction of an appropriate 
counter M'\ Notice that the filter {a^} has a (not necessarily causal) inverse {bj} given 

by 

oo / oo \ — 1 

b jZ i:= $>zM , (7) 

j'=-oo \fe=0 / 

§ Precisely, the atlas containing the single map m : r g R-° — > (r, — r) g R 2 and an orientation 
defined on it. 
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r _1 < \z\ < r. Let C be an upper bound for \x t \ and d be the (minimum) distance 
of the trajectory of {x t } from dM in the maximum norm. For notational convenience 
define x t — yt — for t < 0. Let 

L 

ut ■■= E hvt-j, ( 8 ) 

j=-L 

where L is chosen such that 



\x t - u t \ 



-L-l 



E + E E & 



jOkXt—j—k 



j=L+l I k=0 



< 




C 



< d/2 (9) 

for all integer t. Convergence of the left hand side of Q guarantees that such an L 
exists. {u t } is an approximation of {x t } reconstructed from {y t } using the filter ijSJ). 
Since the approximation error of the time series is at most d/2, so is, in the maximum 
norm, the approximation error of the trajectory. In particular, the topological relation 
between the trajectory and the counter M is not changed by going over from {xt} to 
{u t } (except for some pairs of forward/backward transitions through M, which do not 
contribute to the limit Hence, u>m,u — ^>m,x- 

Now, notice that the ./V-dimensional delay embedding of {ut} can be obtained by 
a linear projection 

(u t , ■ ■ • , ut-{ N ~i)) T = P (yt+L, y t -(L+N-i)) T (10) 

from the 2L + TV-dimcnsional delay embedding of {y t } [HlEIjj with the matrix elements 
of P given by |JSJ. Furthermore, P maps the trajectory of {yt+h} onto the trajectory 
of {ut}. Define the oriented manifold M' such that 

v e M' O Pv 6 M (v g R 2L+N ) (11) 

and dM = PdM' in the obvious way. This guarantees a finite distance of the 
trajectory of {y t } from dM', and there is a one-to-one correspondence between 
transitions of the trajectory of {ut} through M and transitions of the trajectory of 
{yt} through M' . Hence, u>M',y = ^>m,u — ^>M,xi proving Theorem 1. □ 



2. 3. The arbitrariness of the power spectrum, given the topological frequency 

Spectral frequency measures depend on the time series through the power spectral 
density S{u>) alone. By showing that S(u>) is independent of the topological frequency, 
it becomes clear that spectral frequency measures can generally differ arbitrarily from 
the topological frequency, in contrast to what one might intuitively assume (see, e.g., 
reference [H], p 226). 

Theorem 2. Let Sq(l>) > (u) g [— tt, it]) be a symmetric, continuous function, 
< loq < 7r and e > 0. Then there is a time series {yt}, an embedding dimension N 
and a counter M C M. N such that ujM.y = a nd 

| S y (oj) — Sq(uj)\ < e for all u g [— tt, tt], (12) 

where S y (uj) is the power spectral density of {yt}. 

Proof. In order to obtain {y t } as described in Theorem 2, take a time series {x t } 
which oscillates with frequency loq and adjust the spectral density by filtering. A 
suitable time series to start with is given by 

x t — 2 cos(i uiq t + i 4>t) with fa — and <fit+i = <f>t + &eti (13) 
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where {e t } is an equally distributed random sequence of the values —1 and 1. With 
the counter M given by 10, the topological frequency ojm,x = w o is defined when 

0<tf<( W ° . a , !° r °<^ <7r / 2 ' (14) 
[ arcsm(sm wo) for 7r/2 < ujq < n. 

The autocorrelation function of {x t } is (x t x t+T ) t — 2 (cos$) T cos usqT and its spectral 
density 

sin 2 $ fl + cos 2 -!? — 2cos$ COSW COSWn) 

S x {lj) = i 2 ^2 (15) 

7T 1 - e l( ^-"o) cos??| |l - e l ^ +UJ o) cos$\ 

is positive and continuous as required below. In order to see that there is a suitable set 
of filter coefficients {at}, notice that, as an immediate consequence of Theorem 4.4.3 
of reference £Q, there is, for any e > and any two continuous, symmetric spectral 
densities S x (u> ) > OandSo(w) (u> £ [— tt, tt]), a non- negative integer p and a polynomial 
c(z) = 1 + C\z + . . . + c p z p such that 

c{z) ^ for \z\ < 1 (16) 

and, for all u> £ [— tt,tt], 

1 2 5*(w) 

where C = (1 + c 2 + . . . + c 2 )- 1 ^ 1 S(w)/S x (lj) duo. Setting a(z) = C 1 / 2 c(z), 
a = C 1 / 2 , a/,. = C 1 / 2 c / r £ (fc = 1, . . . ,p), and all other = 0, {y t } given by J5J has 
the spectral density S y (ui) — \a(e~ luJ )\ 2 S x (oj) (see, e.g., reference £Q, Theorem 4.4.1) 
and inequality i|17fl implies (|12|l . By (|16|) the filter {afc} satisfies the invertibility 
condition JBJ of Theorem 1. Thus, an appropriate counter M can be obtained 
such that LJM.y — and Theorem 2 is proven. It should be mentioned that when 
Sq(oj) is analytic for real w, Theorem 2 generally holds also with the perfect identity 
S y (ui) = Sq(lj) instead of inequality IT^jl . □ 



C\c(e- tuJ )\ 



maxA o x (\) 



2-4- An Example 

As a demonstration for Theorem 2, consider the time series {yt} shown in figure^. By 
construction, it is a realization of a white-noise process. It was obtained by "bleaching" 
[T2"] a realization of {xt} given by 1)1 3|l with ojci = 1.1 and 9 = 0.4, i.e., the realization 
was filtered such as to transform the known spectral density 1)15(1 into a white spectrum. 
Although all spectral information was lost, ujq can precisely be recovered from {yt}- 
Using an automated search algorithm (to be described elsewhere), a projection matrix 
P (figure^)) is found, such that the projection of the 20D embedding of {y t } into 2D 
yields a trajectory with a nice circular structure and a "hole" in the center (figure^). 
The number of oscillations and the frequency are obvious; ojq is recovered. Any line 
extending from the origin to infinity can serve as a counter M in the 2D projection. 
This can be used to obtain a corresponding counter M' in 20D by a back-projection 
as in jnj. 

A projection with an inadequate P would not yield a different frequency, but 
only a criss-cross kind of trajectory (figure ^i,e), typically with an approximately 
Gaussian distribution of values with a maximum density at the origin. From such a 
representation, no positive topological frequency can be obtained. 

The two-step procedure used here to obtain the counter M' via a counter M in 
2D works for many experimental time series, even though the concept of topological 
frequency is more general. The projector P can then be interpreted as a complex- 
valued filter {fk} with the impulse response function fk = P\,k +iP%k' I n section EPI 
we come back to this point. 
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Figure 1. Demonstration for Theorem 2. (a) A representative segment of the 
white-noise time series yt introduced in section 12.41 (b) the components Pi k 
(solid) and P2 } k (dashed) of a projector from its 20D delay embedding to 2D, and 
(c) a segment of the projected trajectory p(t) of yt, with (u(t),v(t)) T := Pp(t). 
The frequency can be read off. Alternatively (d) another, inappropriate projector, 
and (e) the projected trajectory. 



3. Noisy, weakly nonlinear oscillations 

There are two assumptions upon which Theorem 1 is based - the boundedness of {x t } 
and the finite distance of its trajectory from dM - which are not perfectly satisfied by 
typical noisy processes. Rather, the probability of reaching some point in delay space 
decreases exponentially (or faster) with the distance from some "average" trajectory 
and the inverse noise strength. For many processes the two assumptions and, as a 
consequence, the invariance of the amount under filtering hold therefore only up to an 
exponentially small error. For signals generated by noisy, weakly nonlinear oscillators, 
an analytic estimate of this error shall now be derived. 

Due to the separation of time scales inherent in the weakly nonlinear limit, it is 
more appropriate to work in a continuous-time representation. Consider the noisy, 
weakly nonlinear oscillator described by a complex amplitude A(t) with dynamics 
given by the noisy Landau-Stuart equation J3j 

A = (e + iu )A - (1 + ig t )\A\ 2 A + n(t), (18) 

where e, luq, and gi are real and rj(t) denotes complex, white noise with correlations 

(r,(t)r,(t')) = 0, (r,(t)r,(t'r)=4S(t-t') (19) 
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Figure 2. The power spectrum Sa(^) of A(t) given by equations 1181191 
with e = 2 and gi = 1, obtained from a numerical simulation; and <^ pe ak,A 
(analytic results in reference I15| ~l. compared to the mean frequency a> me an A = 
<^0 - 9i (IAI 2 ) [(|A| 2n ) = 2 n Af- 1 d n Af/de n , see (2g)] and the phase 

frequency oJ p h,A = ^0 — ^ | ^4 1 ^ , defined by equations I2UI21I . and the linear 
frequency u>o- 



[* — complex conjugation, (•) = expectation value]. In a certain sense, this system 
universally describes noisy oscillations in the vicinity of a Hopf bifurcation Jl] . 

3.1. Definitions of frequency 

For the reasons explained above, the topological frequency cannot be defined rigorously 
for noisy, weakly nonlinear oscillators. The customary frequency measures, such as 
the linear frequency cjo, the spectral peak frequency lo vcs ^, the average frequency or 
phase frequency 

Wph.A := (a;*) , Wj := Im {A/ A} , (20) 
and the mean frequency 

_ (loMI 2 ) _ ImjAA*) _ J u;S A (iu)duj 

OWn.A- (|A)a) ^ (|A|2) - J SA{u)duJ W 

will generally (i.e., with gi ^ 0) all yield different values; see figure El 
[Definitions 1|2()I21JI are sometimes restricted to "analytic signals" (Sa (w) = for 
u> < 0) derived from the corresponding real-valued signals Re{^4(£)}. See reference 
[Td] for the history] 

The phase frequency measures the average number of circulations around the 
point A = in phase space per unit time (decompose A(t) = a(t)e 1 ^^ to see this). It 
is a period-counting frequency and the quantity which comes conceptually closest to 
the topological frequency. However, the choice of the point A = can here be justified 
only by symmetry and dynamics [the invariant density pertaining to equation i|18|) has 
an extremum at A = 0], and not by invariance under perturbations. cj poa k and w mean 
are both spectral frequency measures, and the influence of filtering is obvious. But 
how does filtering affect w p h,A? The following considerations lead to a surprisingly 
accurate result. 

3. 2. The effect of filtering on the phase frequency 

The dynamics of A on short time scales St is dominated by the driving noise, and 
the change in A is of the order \SA\ = C(4(5i) 1 / 2 . A band-pass filter of spectral 
width Aw which truncates the tails of the peak corresponding to A in the power 
spectrum suppresses this diffusive motion on time scales Aw -1 , while on longer time 
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scales dynamics change only little. The corresponding deformation of the path of A 
in the complex plane can alter the number of circulations of the origin whenever A 
approaches the origin to less then w (4/Acj) 1 / 2 . At these times |^4| is small and, 
for not too narrow filters, the dynamics of A in its linear range. Thus, the effect of 
broad-band filtering can be estimated by a linear theory. Consider, for a moment, the 
linearized version of equation Ijl8|l . 

i= (e + iu) )A + r){t), (22) 

with r)(t) as above, and assume e < 0. Clearly, w p h,A = ^o- For the phase frequency 
of a complex, Gaussian, linear process B(t) in general, a simple calculation shows 
w p h,B = w m ean,B- This can be used to calculate the phase frequencies of A after 
filtering. Let, for example, B be obtained from A through the primitive band-pass 
filter 



B = {-e x + iwi)B + A, 

which is centered at lo\ with width t\ > 0. Using u> p h,B = w B 
filter theory [2] one obtains 

eiuj -ewi 



(23) 

and elementary 
(24) 



By the argument given above, the shift in phase frequency Slo := Lu p h,B — ^ph,A = 
(ojo — wi) (e/ei) + C(e^ 2 ) is due to the times where \A\ 2 < 4/ei. Since A has a complex 
normal distribution with variance (|^4| 2 ) = — 1/e, this happens about 



P 



\A\ 2 < 



= 1 



exp 



(25) 



of all times. Thus, during these times, the shift in phase frequency is 5w/p[|^4| 2 < 
4/ei] = (lui — luq)/2 + 0(e^ 2 ). Extrapolation to e > and the weakly nonlinear case 
yields 

4 LOl — UJq 



5ui = p 



\ay< 



0{e- 



now with 



where 



1 



10 (el 
exp 



p[\A\ 2 < I 
M = ir 1 ' 2 exp(e 2 /4)[l + erf(e/2)]. 



dl 



h 
N 



(26) 

(27) 
(28) 



Table 1. The shift <5u; num in the phase frequency w p h,A of A(t), obtained 
from simulations of equations 1181191 with gi = 1, after filtering as in 1251 . and 
a comparison with the theoretical estimate 1291 . The data verify 8u ~ e7 , 
~ (lui — loo), and ~ Af~ 1 in this order. 



c 


(^ph,A - <^o) 


(a/i - u ) 








2 


-2.2253 


-2.5 


48 


-0.0117 


-0.0109(16) 


2 


-2.2253 


-2.5 


24 


-0.0235 


-0.0228(16) 


2 


-2.2253 


-2.5 


12 


-0.0469 


-0.0437(15) 


2 


-2.2253 


0.0 


24 


0.0000 


0.0003(17) 


2 


-2.2253 


2.5 


24 


0.0235 


0.0226(16) 





-1.1284 


-2.5 


24 


-0.1175 


-0.1036(47) 


3 


-3.0605 


-2.5 


24 


-0.0063 


-0.0074(13) 


i 


-4.0104 


-2.5 


24 


-0.0011 


-0.0012(11) 
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Equations I|2tjl28fl predict the shift 

5uj = 2e^ 1 Af- 1 (cu 1 -u ) + 0(er 2 ) (29) 

in the phase frequency of A(t) given by i|18H9[l after passing through the filter P5)l. 
A numerical test verifying this result is shown in tabled notice in particular the fast 
decay of Su> as e increases [8uj ~ exp(— 4/e 2 )] and the conditions of Theorem 1 are 
better satisfied. The high accuracy of the result might be understood by observing 
that the crude upper bound \A\ 2 < 4/ei enters the derivation of (|29|) two times, its 
numerical value canceling out. 

3.3. The frequency of real-valued, weakly nonlinear signals 

Experimental signals are real- valued. Assume that, instead of A(t), only a real- valued 
signal x(t) = Kc{A(t) + (higher harmonics)} + (perturbations) is given. The natural 
way to estimate the phase frequency of A(t) then is to construct an approximation 
^(*) = (f * x )(t) °f -Mt) by a convolution of x(t) with a complex-valued filter f(t), 
and to estimate the phase frequency as w p h,A = A- ^he filter f(t) describes the 
combined effect of 2D delay embedding or analytic-signal construction, and filtering to 
eliminate higher harmonics, offsets, aliasing, and external perturbations. The result 
above shows that generally, for w p h,A to be unbiased, the total effect of all these 
transformations should be a complex, symmetric band-pass centered on the linear 
frequency ujq w p h, w P oakO- Otherwise there is a bias which decays as exp(— 4/e 2 ) 
for large e. To the extent that the bias vanishes, the probability of finding values of 
A(t) w A(t) » vanishes, too. Then, a counter M' for x(t) can be obtained along 
the lines of section l2~4l using the filter f{t) - approximated by a time-discrete filter fk 
— for the projection to 2D. Obviously, the corresponding topological frequency u>m',x 
equals c£> p h,A- 

4. Conclusion 

When the spectral density is of genuine interest, forget period counting. But there are 
many real- world applications (e.g., in astronomy jcij, earth science biomedicine 
[SI E] , or engineering 0) where neither the characteristics of the signal pathway 
nor a detailed model of the oscillator are known, and yet a robust measure of the 
frequency or, at least, some robust characterization of the oscillator is sought. Then, 
by Theorem 2, spectral methods miss valuable information. In view of Theorem 1 
and equation concepts such as topological frequency or its little brother, 

phase frequency, are more appropriate. The fractal dimension of the reconstructed 
attractors, an alternative characterization, is typically robust with respect to finite- 
impulse- response filtering only |3EI!> i-e., only if there is a q such that ak = for 
k > q in ©. 

For a practical application of topological frequency, a systematic method to find 
appropriate counters in the typically high-dimensional delay spaces is desirable. Some 
progress in this direction has be made and will be reported elsewhere. 
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